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"^ Abstract 

^^ Multiderivative time integrators have a long history of development, and yet to 

date, they have not been fully explored as a tool for solving partial differential 
I— I equations. In this work, we propose the application of explicit multistage multi- 

■^T derivative time integrators for solving hyperbolic conservation laws. These time 

integrators sit halfway in between popular Runge-Kutta and single-step Taylor 

methods, and therefore, if properly chosen, can leverage advantages from each 

class. Like pure Taylor methods, multiderivative integrators make use of the 

Cd evaluation of higher derivatives of the unknown; like traditional Runge-Kutta 

H methods, multistage multiderivative methods add extra stage values in order 

'~—^ to increase the order of accuracy. We describe a general framework for how 

^_^ these methods can be applied to two separate spatial discretizations: the dis- 

<^ continuous Galerkin method and the finite difference essentially non-oscillatory 

f^ method. Numerical results for a collection of standard test problems are pre- 

^-H sented, and they indicate that the multiderivative time integrators compete well 

0^ with popular strong stability preserving time integrators. 

._J Keywords: Hyperbolic conservation laws, Multiderivative Runge-Kutta, 

f^ Discontinuous Galerkin, Weighted essentially non-oscillatory (WENO), 

C^ High-order Lax-Wendroff. 



1. Introduction 

Effective time integrators form a crucial part of any numerical solver for 
partial differential equations (PDE). In this work, we revisit classical ordinary 
differential equation (ODE) solvers known as multiderivative (Obreshkoff) meth- 
ods. These methods include evaluating higher derivatives of the function in their 
framework, which reduces the total number of stages that would be required for 
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achieving the desired order of accuracy. The motivation for this work comes 
from a desire to develop numerical methods that are best suited for modern 
computer architecture, where recent trends demand considering floating point 
operations (FLOPs) as inexpensive, and memory access and storage as a scarce 
resource that should be reduced at all costs. Our current focus is to demon- 
strate how multiderivative Runge-Kutta methods can be used as an effective 
time integrator for solving hyperbolic conservation laws. 

A conservation law is a PDE defined by a flux function R of the form 

9.t + Vx • i?(g) = 0, g(0,x) = go(x), x e 17 C M'*. (1) 

In dimension d, this problem is deflned by m equations with initial conditions 
prescribed qq : M.'^ — > M™. The unknown quantity q is a function of time as well 
as space: g : E+ x M"^ ^ M™, and the flux function R : M'" ^ M'' x M'" defines a 
physical quantity of interest, with each component from physical space denoted 
as R = [/'■^^ f^^\ ■ ■ ■ , Z'-''-']"^. The system is hyperbolic if the matrix 

n^ ' -Q— (q) + ri^ ' -g— (q) + --- n^ ' -Q— (q) 

is diagonalizable for every n e M'' satisfying ||n|| = 1 and q in the domain of 
interest. 

Numerical methods for solving (fTj) require a discretization technique for space 
as well as a (possibly coupled) discretization technique for time. The vast ma- 
jority of time stepping discretizations fall into one of two distinct categories: 

• Method of lines formulation. 

A method of lines (MOL) solver for (IT]) separates the spatial discretization 
from the time evolution. Starting with q,t = — Vx ■ R^q), one defines 
a spatial discretization q'^ of the continuous variable q, which could be 
tracking point values (finite difference, spectral methods), cell averages 
(finite volume methods), or coefficients of basis functions (finite element 
methods). This operation defines a function C{q^) = — Vx • R{q^), which 
in turn defines a large ODE system of the form: 

q':t - C{q% (2) 

Once this discretization has been parsed, one may apply an appropriate 
ODE integrator to this problem: for hyperbolic conservation laws, explicit 
time-stepping methods are usually preferred. 

• Lax-Wendroff (Taylor) formulation. 

The Lax-Wendroff [1 procedure is a numerical scheme that updates the 
solution using finitely many terms from the Taylor series of the function. 
Here, one first expands q(t,x) in time about t = t": 

(t -f^V 
q(t, x) = q" + {t- f^)q^ + ^-^^q'^t + ' ' ' , (3) 



and then each time derivative is replaced with a spatial derivative via the 
Cauchy-Kowalewski procedure: 

9.t--Vx-i?(g), (4) 

q,tt = -Vx • R{q)t = -Vx • {R'{q)q^t) = Vx • [R' {q) ■ Vx • R{q)) , 



Further derivatives are required for higher order variants, and of course, 
one still needs to choose a spatial discretization. Inserting t = t" + Ai 
produces a single-stage, single-step method. In addition, this is funda- 
mentally different than the MOL discretization, because the physical and 
temporal variables are intimately intertwined through the choice of the 
spatial discretization of the right hand side of ((sl) . 

The Lax-Wendroff procedure is much older than either Lax or Wendroff. 
Indeed, a more appropriate name would be the Cauchy-Kowalewski procedure. 
Augustin-Louis Cauchy (1789 - 1857) and Gerhard Kowalewski (1876 - 1950) 
each sought methods that could aid them in proving existence and uniqueness for 
solutions to PDEs. Their combined method, known as the Cauchy-Kowalewski 
procedure, was presented in equation (Wl), was derived from Brook Taylor's 
method, who invented equation (Is]) in the 1700's. For a modern (mid-20"' 
century) proof and review of the Cauchy-Kowalewski procedure see Friedman 
[2 or Fusaro [3] and references therein. In 1960, Peter Lax and Burton Wen- 
droff [1] realized the Cauchy-Kowalewski procedure could be used as a numerical 
method. They started with the theoretical groundwork developed by Cauchy 
and Kowalewski and derived a numerical scheme for solving PDEs. Therefore, 
this entire procedure is often cited as the Lax-Wendroff method within the nu- 
merical analysis community. 

The original Lax-Wendroff method was a second-order numerical discretiza- 
tion of the Cauchy-Kowalewski procedure, and over the past decade, much work 
has been put forth to define high-order variants of this method. Perhaps the 
most direct extension of the Lax-Wendroff method is derived by looking at suc- 
cessively higher order modified equations, and adding extra terms to account 
for further derivatives. Daru and Tenaud [3| explored high-order monotonicity 
preserving single-step methods, and they derived TVD flux limiters to control 
spurious oscillations. In 2003, Jianxian Qiu and his collaborators demonstrated 
a high-order extension of the Lax-Wendroff procedure using finite difference 
weighted essentially non-oscillatory (WENO) methods |S]. Later on, they ap- 
plied the same procedure to Hamilton- Jacobi systems as well as the shallow 
water equations [6l [7] . Additionally, Qiu, Dunibser and Shu showed how to ap- 
ply the Lax-Wendroff scheme to the discontinuous Galerkin (DC) method [8], 
and shortly thereafter, Dumbser and Munz followed a similar procedure for con- 
structing DG methods to arbitrarily high orders of accuracy using generalized 
Riemann solvers ^ . High-order versions of a Lax-Wendroff type of discontin- 
uous Galerkin scheme have been investigated for ideal magnetohydrodynamic 
equations |10j . and explorations into various numerical flux functions for the 



Lax-WendrofF DG method has also been carried out [TT]. It has already been 
noted that high-order schemes with Lax-WendrofF type time discretizations can 
be implemented to carry a low-memory footprint |T2 . 

Much of the difficulty when constructing high-order versions of the Lax- 
WendrofF scheme comes from the necessity of defining higher spatial derivatives 
of the solution. After producing the Jacobian of the flux Function, the next time 
derivative produces the Hessian of the flux Function. Further derivatives require 
tensors which grow vastly in size, and thereFore, one of the primary concerns 
with a high-order Lax-WendrofF method is the burden of implementing higher 
derivatives, especially in higher dimensions. We believe that this is the primary 
reason we haven't seen high-order Lax-Wendroff methods gain traction in large 
scale codes: many developers and users oF these codes don't have access to these 
derivatives, or rather, they do not have the resources to develop the proper 
interFaces to build them into their codes. The trade-ofF For computing higher 
derivatives is that one ends up with a single-stage, single-step method, which 
can be implemented with the optimal, low-memory footprint for hyperbolic 
problems. 

On the other extreme lies a pure MOL Formulation of the PDE, which given 
its portability, is the predominant method in the community. This Formulation 
is incredibly attractive From a soFtware engineering perspective: one can en- 
vision developing a code that completely decouples the ODE technology From 
the spatial discretization, whereas Taylor methods require a recognition oF the 
particular spatial discretization. In addition, the MOL formulation invites de- 
velopers to concentrate efForts on ODE integrators as a separate entity From 
the PDE. However, this idealization is not possible given that the scheme is 
intended to solve a PDE, and thereFore one needs to respect the choice oF spa- 
tial discretization not only when selecting an ODE integrator, but also when 
developing one. 

In this work, we advocate the use of multistage multiderivative integrators For 
solving hyperbolic conservation laws. We believe common ground can be Found 
to leverage advantages From MOL Formulations and pure Taylor methods. We 
argue that a wise choice oF the correct multiderivative integrator can accomplish 
the Following: 

• High-order accuracy [3^'^ -order or higher]. Explicit multiderivative 
schemes can be constructed to arbitrarily high orders of accuracy. We 
focus on a single fourth-order method as our demonstrative example. 

• Portability. Hyperbolic codes usually have access to the eigen-decomposition 
oF their problem. Multistage multiderivative integrators that stop at the 
second derivative don't require For anything above and beyond this, and 
thereFore, they are more portable than pure Lax-Wendroff methods. 

• Low-storage. Multiderivative integrators carry a small memory foot- 
print. By design, these integrators exchange storage For extra FLOPs 
in order to attain high-order accuracy. This Feature is a desirable trait 



for high performance computing given than the current trend is towards 
inexpensive FLOPs and expensive memory. 

The primary purpose of the present work is to demonstrate how multistage 
multiderivative integrators can be used to solve PDEs. Given the plethora of 
combinations concerning the number of stages and number of derivatives, our 
current work presents a single model example. This example is general enough 
to serve as the foundation for an investigation into any of the multitude of 
combinations involving extra stages and extra derivatives. 

The outline of this paper is as follows: we begin in S|2]with a historical review 
of multiderivative integrators. In S|3J we describe the finite difference WENO 
scheme, and we demonstrate how multiderivative technology can be applied to 
the WENO framework. In Q we continue by looking at multiderivative integra- 
tors for the discontinuous Galerkin method. In SjSlwe present numerical results 
for our numerous numerical test problems, and in fp] we draw up conclusions 
and point to future work. 

2. High-order explicit multiderivative ODE integrators 

Multistage multiderivative integrators for PDEs require a blend of both the 
the method of lines (MOL) formulation as well as the Lax-Wendroff formulation 
of (IT]), and in addition, one needs to select a method for the spatial discretiza- 
tion. We begin our description of multistage multiderivative PDE solvers with 
a historical overview of these methods within the confines of ODEs. In particu- 
lar, we focus on explicit multiderivative Runge-Kutta integrators, which include 
single-derivative methods (e.g. classical Runge-Kutta) as well as single-stage 
methods (e.g. Taylor) as special cases. We begin with a review of multideriva- 
tive methods in |2.1[ and then continue in j j2.2| by defining a large class of 
explicit multiderivative Runge-Kutta methods. In |2.3[ we present the model 
example used for the entirety of this work. 

2.1. Multiderivative methods: a review 

Numerical methods using multiderivative technology have a long history dat- 
ing back to at least the early 1940's, and some of the pioneering work for explicit 
schemes share a common ancestry with implicit schemes. In 1963, Stroud and 
Stancu [T3] applied the quadrature method of Turan [T3] , which generated an 
implicit, high-order multiderivative ODE solver. Prior to Turan's work, in 1940, 
Obreshkoff ^15] derived discrete quadrature formulae for integrating functions, 
and much like Turan did, Obreshkoff used extra derivatives of the function for 
his quadrature rules. When extra derivatives are included, one can obtain meth- 
ods with remarkable properties for the numerical integration of ODEs, including 
high-order accuracy involving fewer quadrature points than would otherwise be 
required. In 1972, Kastlunger and Wanner [TB] used the theory of Butcher 
trees to show that Turan's quadrature method could be written as an implicit 
multiderivative Runge-Kutta scheme. The following year, Hairer and Wanner 
[17] defined "multistep multistage multiderivative methods" , which to date, has 



stood the test of time as being an incredibly broad categorical definition of nu- 
merical methods for solving ODEs. In particular, their definition contains all 
Runge-Kutta and all linear multistep methods as well as additional combina- 
tions, including so-called general linear methods (c.f. John Butcher's extensive 
review papers [TH [121 HO] for a description of general linear methods). The 
textbooks of Hairer, N0rsett, and Wanner [5T1[25] contain excellent references, 
and should be required reading for anyone interested at looking into the history 
of multiderivative methods. 

Our current focus is on explicit versions of multiderivative Runge-Kutta 
schemes, which needless to say, also have a long history of development. De- 
spite their age, these methods have seen little to no attention outside the ODE 
community, yet given the direction of modern computer architecture, many 
of these methods may see use for solving PDEs in the near future. In 1952, 
Rudolf Zurmiihl j23j investigated multiderivative Runge-Kutta integrators, and 
later on, Erwin Fehlberg [211 US] derived an explicit multiderivative Runge- 
Kutta methods. Early versions of Fehlberg's method applied a single-derivative 
Runge-Kutta method to the modified variable that is constructed by subtract- 
ing out m-derivatives of the original variable. A decade later, Kastlunger and 
Wanner [26] extended Butcher's method to multiderivative Runge-Kutta meth- 
ods. In their work, they defined the order conditions for the coefficients in a 
multiderivative Butcher tableau, and in addition, they showed that Fehlberg's 
method [211 123 <^^^ be written as a multiderivative Runge-Kutta process with 
TO-derivatives taken at a single node. During that same decade, Shintani [2711^ 
worked on multiderivative Runge-Kutta methods. Also in the 1970's, Bettis and 
Horn [29] revisited Fehlberg's scheme and reformulated it as an embedded Taylor 
method: for their celestial mechanics problem, they describe how the necessary 
Taylor series coefficients can be generated with little to no additional cost. A 
decade later and unaware of the full history of the methods, Mutsui [SU] also 
worked on Runge-Kutta methods that leveraged extra information with extra 
derivatives. 

The most recent work on explicit multiderivative integrators appears to focus 
on redefining order conditions and demonstrating examples of methods from this 
class, much of which has been carried out independently from previous work. 
In 1986, Gekeler and Widmann [3T] used the theory of Butcher trees to define 
the correct order conditions for multiderivative Runge-Kutta methods. In their 
work, they presented families of methods with orders ranging between four and 
seven. Goeken and Johnson [321 IS3] were unaware of the long standing history 
of explicit multiderivative methods, and they derived their own versions of ex- 
plicit methods that are sub-optimal. Within the past decade, Yoshida and Ono 
[34] [35] and Chan and Tsai [36] used the theory of Butcher trees to define order 
conditions and presented numerous examples. The primary difference between 
the latter two works is the following: Chan and Tsai used multiple stages for 
their methods, and they restrict their attention to using two derivatives of the 
unknown function; Yoshida and Ono restrict their attention to two-stage meth- 
ods, and they admit arbitrarily many derivatives of the unknown function to 
be evaluated at every quadrature point. In other very recent work, Nguyen-Ba, 



Bozic, Kengne and Vaillancourt [37] derived a nine-stage explicit multiderivative 
Runge-Kutta scheme that makes use of extra derivatives at a single quadrature 
point only, much like the schemes Fehlberg originally investigated [Ml ESj • In 
doing so, the order conditions become simpler to navigate. 

For the purposes of solving hyperbolic conservation laws, we view using at 
most two-derivatives of the function as the most appropriate choice given the 
opportunity to retain portable code. Investigations into methods using extra 
derivatives would make for interesting future research. 

2.8. Multistage multiderivative methods: some definitions 
Consider a system of ODEs, defined by 

y^Liy), y(0) = yo, t > 0. (5) 

We use the letter L in place of / to avoid conflict with the flux function defined 



later on in equation (13). Without loss of generality, we assume the system 
is autonomous. Multiderivative methods make use of extra derivatives of (Is]), 
which obey the recursive relationship 

L^'^Hy)^^ — y = % — Hy), mez>i. (6) 



(m+1) 



dy " dy 

Definition 2.1. A multiderivative Runge-Kutta scheme (MDRK) with s-stages 
and r-derivatives is any update of the form 

yn+i = 2/« + E E 6!™^Af"L(™-i)(2/(')), (7) 

2^1 m— 1 

where intermediate stage values are given by 

y'-'^ =yn + i2i2 a'i^^^At"^L'-"^-'Hy(^y). (S) 

j — 1 771 — 1 

//aj" — whenever i < j , the method is explicit, otherwise it is implicit. 
We remark that both Taylor and traditional Runge-Kutta methods are spe- 



cial cases Definition 2.1 setting r = 1 produces traditional Runge-Kutta meth- 
ods, and setting s = 1 produces the Taylor class of methods. 

Our definition is an equivalent, yet distinctly different version of what can 
be found in other sources (c.f. [H]). It is possible to define interm ediat e stages 
through defining and saving L('"^(y(*^), but we prefer Definition 2.1 because 



of the potential for a low storage implementation, at the cost of recomputing 
previously observed values. We claim that this formulation is precisely what 
makes multiderivative integrators competitive with current state of the art ODE 
integrators, under the assumption that FLOPs are inexpensive, memory usage 
and communication is to be reduced at all costs, and that the relative size of y 
versus /(y) and L'^"^\y) are all on par with each other. 

For hyperbolic conservation laws, we consider methods that use at most 
two-derivatives to be the most portable. 



Definition 2.2. An s -stage, two- derivative Runge-Kutta (TDRK) scheme is 
any update of the form 



y^+,=y,, + AtJ2b^^Hy^'^) + ^t'J2^^'^^{i'''^. (9) 



1=1 i=l 

where interm,ediate stage values are given by 






If a"^ — for all i < j, the method is explicit. 



«j 



Before presenting an example of a method from this class, we would like to 
draw some comparisons between popular time integrators and the multideriva- 
tive Runge-Kutta methods. Our aim is to discuss advantages each method has 
for being coupled with numerical PDE solvers, and in particular, we would like 
to focus on which methods work well with new computer architectures. 

Traditional Runge-Kutta methods are far and wide the most popular for 
solving hyperbolic conservation laws, yet we see room for improvement given 
the current direction of computer architecture. Runge-Kutta methods are easy 
to implement, and therefore, they are the most portable of all of the methods 
discussed here. They are self-starting and can easily change their time step 
size, which is an important characteristic to have for solving hyperbolic conser- 
vation laws. In addition, when compared with their natural counterpart, the 
Adams family of methods (e.g. linear multistep methods), Runge-Kutta meth- 
ods have stability regions that are more favorable for hyperbolic problems. For 
example, on a purely oscillatory problem, the maximum stable time step for 
classical fourth-order Runge-Kutta is given by \z\ < a/S « 2.8, where z = AAi 
is purely imaginary. For the same cost and identical storage, one would be able 
to take four time steps with fourth-order Adams Bashforth. Even after rescal- 
ing, the maximum stable time step for the equivalent Adams method would be 
restricted to |z| ^ 1.72. It would seem that Runge-Kutta methods are ideally 
suited for solving hyperbolic conservation laws. They can be derived to require 
low-storage |38l ESI HQl IS] , can be designed to acquire strong stability preserv- 
ing (SSP) properties [1111131133], and are very portable, especially given that 
they are self starting. However, traditional Runge-Kutta methods are not opti- 
mal with their memory usage, and to date, even the low-storage Runge-Kutta 
methods require many stages, and therefore they may require considerable com- 
munication overhead when compared to pure Taylor schemes. 

Taylor methods lie on the other extreme: we claim that they can be imple- 
mented to have optimally low-storage for hyperbolic problems, and can contain 
minimal communication overhead. However, pure Taylor methods are the least 
portable of the time integrators discussed here. In order to implement a high- 
order Taylor (e.g. Lax-Wendroff) method for solving a PDE, one needs to have 
access to high derivatives of the unknown, which puts them out of reach from 



many scientists. We recognize that this can certainly be done for very compH- 
cated problems (TU], but it's difficult to convince users of legacy codes to modify 
them in order to reach high-order accuracy. On the plus side, Taylor methods 
naturally contain favorable stability regions for hyperbolic conservation laws, 
and given that they're single-step methods, they have nominal communication 
overhead. 

We view multiderivative Runge-Kutta methods as sitting halfway in between 
traditional Runge-Kutta and pure Taylor methods, and therefore, can retain 
desirable qualities within each method. In order to retain portability, we view 
multiderivative Runge-Kutta methods that use at most two-derivatives are op- 
timal for hyperbolic conservation laws, especially given that most codes already 
have access to, or at least users would be willing to implement the Jacobian 
of the flux function. Beyond two-derivatives, we argue the "many" -derivative 
Runge-Kutta methods start to lose their portability. 

2.3. Multistage multiderivative methods: an exam,ple 

The example used for the remainder of this paper will now be presented. 
There is a unique combination for an s = 2-stage method that produces fourth- 
order accuracy [3S1 [351 IMj • This method, which we refer to as TDRK4, can be 

written as 





At ^ {At 


y* 

Vn+l 


Vn 1 r, ^n 1 

= Un + AiL„ + — 



2 
3 



-9{Vn). (11) 



-{9{yn) + My*)) 



where we have renamed the first time derivative of L with g{y) — ^L{y). 

Our goal is to describe how to implement the large class of explicit multistage 
multiderivative methods for solving PDEs, and therefore, for simplicity of ex- 
position, we use this method as our canonical example of an explicit multistage 
multiderivative integrator. This method has a stability region identical to clas- 
sical fourth-order Runge-Kutta, and therefore, it will have favorable properties 
for solving hyperbolic conservation laws. In addition, we claim that a complete 
description for this scheme will provide the necessary mechanisms for extension 
and investigation into integrators containing extra stages. These integrators can 
be constructed to be even higher order accurate and contain favorable stability 
regions. For example, all of the methods presented in Chan and Tsai [3^ can 
be implemented using our description with a straight-forward extension of what 
follows. 

The notation used for the remainder of this paper will focus on how to define 
a two-stage, two-derivative method that computes the following quantity: 

y = Vn + {aAtLn + PAt^giVn)) + {a*AtL* + (3*At^g{y*)) . (12) 

In this equation, y could be the full update, as in y = Vn+i from equation (|9| 
or a stage value y = y^*-' from equation (10). Observe that the TDRK4 method 



in equation (11) can be cast in the form of equation (12 1, with the first stage 
given by 

a = 1/2, P = 1/8, «*=/?*= 0, 

and the update given by 

a = 1, /3 = 1/6, a* = and /3* = 1/3. 



We prefer introducing a and /3 over a^™ and 6^ from Definitions 



2.1 



and 



2.2 



because these letters delete unnecessary indices and the upcoming descriptions 
for the PDE methods will introduce further indices that can become cumber- 



some. Equation (12) restricts our attention to two-derivative methods that rely 
on the previous two stage values only (e.g. low-storage), but we remark that ex- 
tensions to multi-stage, 'many'-derivative methods follow by adding extra terms 



to equation (12) 



This completes our description of the multiderivative Runge-Kutta scheme. 



and it bears repeating that without loss of generality, we will use method (11) 
as our canonical example of a method from this class. Given that the focus of 
this work is how to implement this scheme for solving hyperbolic PDEs, we still 
need to describe how to discretize in space. 

3. The finite difference WENO method 

The finite difference weighted essentially non-oscillatory (WENO) method 
has many variations and a long history of development. The original method was 
developed by Shu and Osher [i?! [33] , and later analyzed and further developed 
by Shu and his collaborators [iSl H?! |3H] ■ For a recent comprehensive review of 
the many variations of WENO schemes, see Chi- Wang Shu's extensive review 
paper [49] . In this work, we consider the fifth-order WENO-Z scheme [501 [HI 
[52] , that is an improvement of the classical fifth-order Jiang and Shu (WENO- 
JS) scheme [1^. The underlying choice of the reconstruction procedure is not 
central to this work and our description will be generic enough to accommodate 
most variations of the classical WENO method. However, given that there are 
many options, we explain the minimal details necessary to reproduce the present 
work. 

3.1. The finite difference WENO method: MOL formulation 

We begin our description with a reduction of (II]) to a ID conservation law: 

g,t + /(g). = 0, q{Q,x)^qa{x), x£n^[a,b]. (13) 

Here, we follow the standard convention of naming our flux function / in place of 



R. Moreover, (13) is hyperbolic if the Jacobian f'{q) = RKR ^ is diagonalizable 



with real eigenvalues for all q in the domain of interest. 
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The spatial discretization for a finite difference method seeks a point-wise 



approximation to the exact solution of ( 13 ) at a finite collection of points. We 



start with a uniform discretization of [a, b] using rux points: 

Ax = {b- a)/m^, Xi = a + {i - l/2)Ax, i e {1,2, . . . ,m^} , 

and we seek values qi that approximate the exact solution at each grid point, 
qi{t) « q{t,Xi). 



In order to write (13) in discrete flux-difference form so we can have a con- 
servative methocrj we begin by defining an implicit sliding function h through 

^ l.x+\x/2 

f{q{t,x))^— h{t,x)dx. (14) 

^X Jx-Ax/2 

With this definition in place, we have the nice result that 

— {q{t, Xi)) = -^[h {t, 2:i+i/2) - h {t, a;,,_i/2)] , (15) 

which is precisely what is needed to define a discrete flux difference formulation 
of gt = —fx- The MOL formulation for the finite difference scheme defines 
interpolated values ^^+1/2 that approximate h (t, Xi^i/2) to high-order accuracy. 



1 /9f 

— (/.,+i/2 - h.-i/2) = ^ {q{t, x,)) + {Ax'') 



Moreover, an astute observation means that h need never be computed |171 HU] '■ 
define cell averages of h[t,x) as 

hi := — / h{t,x)dx, 

then observe that f{qi{t)) — hi, and therefore point values of / can be inter- 
preted as cell averages of h. After defining an appropriate interpolating algo- 
rithm for producing high-order interface values /ii+1/2 from cell averages hi, the 
full MOL formulation is given by, 

J^l^i*) = -^ ('*«+i/2 - ^^-1/2) ■ (16) 

This scheme is automatically conservative as it is written in flux difference form. 



Usually one applies a high-order explicit Runge-Kutta integrator to (16), 
which results in what is normally called the Runge-Kutta WENO (RK-WENO) 
method. 



^A finite difference metiiod is conservative if the metfiod satisfies -jh (X]i Qii^)) = on a 
periodic (or infinite) domain. 



11 



3.1.1. The finite difference WENO method: the reconstruction procedure 

A conservative reconstruction procedure requires a single value /ii+1/2 for 



equation (16) to be defined at each grid interface. In the ensuing discussion, 
we suppress the time dependence of h, and assume that we have known cell 
averages hi for a function h = h{x). 

The fifth-order WENO method uses a five point stencil shifted to the left or 
right of the interface: 

'*j+i/2 -^ WENOb~\hi^i,hi,hi+i,hi+2,K+3]- 

Here, we define coefficients for the function WENO'b'^ , and note that by symme- 
try, the reconstruction procedure weighted in the other direction can be observed 
by hipping the stencil: 

WENOb-[h-iX. ■ ■ ■ h+3] ■■= WEN05+[h,+3, h+2, ■ ■ ■ h-i]- 

Three sub-stcncils Sq = {hi^2,hi-.i,hi}, Si — {hi-i,hi,hi^i}, and ^2 = 
{hi, /li+i, /iz-i-2} uniquely define three quadratic polynomials Pj{x) that have the 
same cell averages for each element in their stencil. Each polynomial defines a 
competing value for h{xi+i/2) with /i-+\/2 = Pji^i+i/2)- 

\>i/2 = 3^^-2 - g'l^-i + -^h,, (17) 

/i.l+1/2 = -g^i-i + g^i +3/i*+i, (18) 

hf+i/2^ ^h + -^h+i~-^h+2- (19) 

The linear weights jj — {1/10, 3/5, 3/10} are defined as the unique linear 



combination of equations ( 17)-([19| that yields a fifth-order accurate point value 
for h{xi+i/2): 

^«+i/2 = 7o/j1+\/2 + 7l/i!+i/2 + 72/l!+i/2- (20) 

The WENO procedure replaces the linear weights 7^ with nonlinear weights ujj 
that are necessary in regions with strong shocks. 

The Jiang and Shu smoothness indicators Pj place a quantitative measure 
on the smoothness of each stencil based on a Sobolev norm: 

/3,:=5:Ax--W N^P.(-) dx, 

which for our fifth order method, are given by 

/3o = — {h,_2 - 2/i»-i + K)^ + - {h,_2 - -ih-i + ihf , 

13- --2I- -2 

A = ^ {h^^l - 2hi + h,+i) + - (/i,_i - h,+i) , (21) 

13- - -2I- - -2 

/?2 = T- {hi - 2h.,+i + hi+2) + T (3/ii - 4/ii+i -t- hi+2) ■ 
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The non-linear WENO-Z weights u!^ are a slight modification of the classical 
weights LUj. The new weights require the computation of a single additional 
parameter r^ = |/32 — /?oh 

We use the power parameter p = 2 and regularization parameter e — 10^^^ for 
all of our simulations. 

With these definitions in place, the final interpolated value is defined as 

WEN05+[h_2, ■ ■ ■ k+2] ■■= ^oh%2 + ^i^?i/2 + ^2/»-+\/2- 

3.1.2. The finite difference WENO method: MOL formulation for systems 

For a system of conservations, the reconstruction procedure described in 
§3.1.1| is carried out locally on each of the scalar characteristic variables: 

Finite difTerence WENO procedure for ID systems 

1. For each i, evaluate fi = f{qi) at each mesh point and compute average 
values of q at the half grid points 

<ll-l/2^1^{<l^+<l^-l)■ (22) 



Roe averages [53] may be used in place of ( 22 ), but all the results presented 
in this work use this arithmetic average. 

2. Compute the left and right eigenvalue decomposition of f'{q) = RAR^^ 
at the half-grid points: 

Ri-1/2 = R{q*-i/2)^ ^j^i/2 = ^~ (9*-i/2)- (23) 

Compute a := max^ |/'('?i*_i/2)l &s the fastest wave speed in the entire 
system. For stability, we follow the common practice of increasing a by 
exactly 10% in order to completely contain the fastest wave speed, that 
is, we set a = 1.1 • max^ |/'('7*_i/2)l- The value |/'(9)| is the maximum 
absolute value of all eigenvalues of the Jacobian, f'{q). 

3. For each i, determine the weighted ENO stencil {i + r} surrounding i. In 
fifth order WENO, the fuU stencil is given by r e {-3,-2,-1,0,1,2}. 
Project each qi^j. and flux values /i+r onto the characteristic variables 
using R~\/^: 

Wi+r = R~\/2 ■ li+ri (24) 

gi+r = Ri-i/2 ' fi+r- (25) 

Apply Lax-Friedrichs flux splitting on gt+r'- 
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As an alternative, one could use a local wave speed. Note that for each 
component, these definitions automatically satisfy -^ > and -^ < 0. 

4. Perform a WENO reconstruction on the characteristic variables. Use the 
stencil which uses an extra point on the upwind direction for defining g^ : 

.9+1/2 - WEN05+ [g+ 3, .g+ 2, .9+ i, 5+, .9++i] , (27) 

97-1/2 = WENO5- [9^12^9^-1,9^,9^+1,9^+2] ■ (28) 

Define 3i_i/2 := ff+ 1/2 + 5,"-i/2- 

5. Using same projection matrix, Ri-1/2, project characteristic variables 
back onto the conserved variables: 

ft-i/2 := Ri-1/2 ■ 9i-i/2- (29) 

3.2. The finite difference WENO method: a multiderivative formulation 

In order to implement multiderivative methods into the WENO framework, 
we closely follow previous work on Lax-Wendroff WENO methods [5], but in 
place of relying on a single-step Taylor series, we we build intermediate stages 



of the form given by (12 1. Consider two 'stages', g" and q* that provide a 
pointwise approximation to the exact solution. If we apply Q to each g" and 



q* and insert the result into (12), we see that the starting point for putting 



together a multiderivative WENO integrator is: 

q = 9"~«At/(g"),-a*Ai/(g*), (30) 

- 13 At' {fiq'W^, - f3*^t' {ni*Kt),, ■ 

Note that we have retained the value q^t in place of substituting q^t = —f.x 
into the last two terms. The complete multiderivative procedure is given by the 
following procedure: 

Multistage multiderivative WENO procedure 

1. Given two pointwise approximations g" and g*, perform a single WENO 
reconstruction on each piece to construct the following values: 

It.t -^ ~-^ yfi+i/2~ fi-i/2) ^ 1i,t -^ ^-^ yfi+1/2^ fi-1/2) ■ (31) 

2. Define the pointwise values Gf :— /' (gf ) g"^ and G* :— f (g,*) g,*^. Note 
that we do not need to decompose g onto the characteristic variables. 

3. For each G" and G* in Step 2, compute a finite difference approximation 
to G"!^ and G%. Given that this term inherits an extra factor of At, we 
can use the fourth-order centered finite difference: 

D'xGi :— — — — (Gi-2 ~ 8Gi_i + 8Gi+i — Gi+2) ■ (32) 

12/ax 
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Using a centered stencil means that the method will automatically be 



conservative. The complete WENO discretization of (30) is now, 

q, =q^ + aAtqlt + a* Atql^ - pAt^D^G'l - p*At^D,G*. (33) 

We note that further stages can be accommodated by simply adding more 
terms into ( |33[ ). In addition, further derivatives can be included, but one would 
need to revisit Q before inserting these terms. Higher derivatives can be ap- 
proximated using smaller centered finite difference stencils because they get 
multiplied by increasing powers of At (c.f. [5] for further details). Equation 



(33) finished the spatial discretization of il2n, which can be applied twice to 



construct the unique, two-stage, two-derivative Runge-Kutta method, TDRK4 



in (111 



Qt ^<lt 



(/r+i/2 - /r-1/2) - ^^^ {D^G.,) ■ (34a) 



2Ax 

At ( ;„ ,„ x Ai2 



9?+' = ^r - ^ (/r+1/2 - /r-1/2) - -^ {D.G, + 2D,G:) . (34b) 

Numerous examples of the proposed WENO scheme are provided in Sj5J We 
compare the proposed time integrator with classical fourth-order Runge-Kutta 
(RK4), as well as the third-order strong stability preserving (SSP) method of 
Shu and Osher [44'. Before presenting results, we first describe an implementa- 
tion of a multistage, multiderivative discontinuous Galerkin method. 

4. The discontinuous Galerkin (DG) method 

The discontinuous-Galerkin (DG) method dates back to 1973 when Reed & 
Hill [54] developed the scheme for solving a neutron transport equation. The 
theoretical framework for the DG method was solidified by Bernardo Cockburn 
and Chi- Wang Shu through a lengthy series of papers. We refer the reader to 
their extensive review article and references therein [55]. In this section, we 
define the notation used for the remainder of this paper and provide minimal 
details necessary for reproducing this body of work. This section focuses on 
describing the Runge-Kutta discontinuous Galerkin (RKDG) scheme, and in 



^4.2 we introduce the multiderivative technology to the DG framework. We 
use similar notation to that was previously introduced [56J, and for further 
details on the DG method, we direct the reader the references (e.g. [55} [57]). 

4-.1. The discontinuous Galerkin (DG) method: a MOL formulation 

The DG method solves a discretization of the weak formulation of the hyper- 



bolic conservation law ( 13 1. The continuous weak formulation can be realized by 
multiplying (13 1 with a test function Lp, and integrating by parts over a control 
volume T = [xi, Xr]'- 

ipq^t dx= ip^^ f(q) dx - (^ipfiq) \^^ - (fifiq) |^ J . (35) 
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We begin our discretization by defining a grid containing nix cells for the do- 
main [a, b], each of whose width is Ax = {b — a)/mx- The i*'' grid cell is denoted 
byTi = [a;i_i/2,.Ti+i/2], where the cell edges are given by a:i_i/2 = a + {i — l)Aa:;, 
for i = 1, 2. . . . , rrix + 1, and the cell centers are given by Xi — a + {i — l/2)Aa;, 
for i = 1,2,..., rrix- For simplicity of exposition, we will restrict our attention 
to a uniform grid. On this grid we define the broken finite element space 

W'' = {w'' e L°^in) : w'^lre PP, vr e %} , (36) 

where h — Ax. The above expression means that on each element T, w^ will be 
a polynomial of degree at most p, and no continuity is assumed across element 
edges. Each element Ti can be mapped to the canonical element ^ G [^l, 1] via 
the linear transformation 

x = x, + ^^. (37) 

Note that after a change of variables, spatial derivatives obey the following rule: 
-s- — -K-4c- For the canonical element, we construct a set of basis functions 

ax i\x at, ' 

that are orthonormal with respect to the following inner product: 



(^(^^W):^i|\w(O^W(Ode 



"J Ik, 



where 5ik is the Kronecker delta function. This defines the Legendre basis 
functions: 



We consider approximate solutions of the hyperbolic conservation law ( 13 1 
that are defined by a finite set of coefficients \ Q\ \ of the basis functions 
{(^^'^•'j. When restricted to a single cell 7^, the approximate solution q'^ is 



M 



q\t,x) := gf (i,e) = ^ QW(i) ^W(e), (38) 

k=l 

where M is the order of accuracy. The initial conditions are determined from 
the i^-projection of q*^{{),x) onto the basis functions, 

Qf\Q):=(q1{Q,0,V^^\i)), (39) 

which we evaluate using M standard Gaussian quadrature points. 



The semi-discrete weak formulation of ( 13 ) is given by discretizing ( 35 ) with 



a finite set of basis functions and a finite set of control volumes. Inserting our 



discrete spatial representation ( 38 ) into the continuous weak formulation ( 35 1 
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and using (p — ip^^' for our test function, we produce the semi-discrete weak 
formulation 



(g5(t,0,^('=HO) = ^(/(9f(t,0),^f'(e)) (40) 



where the flux values f^{q^(t,Xi-i/2)) will be defined through an appropriate 
Riemann solver. These terms will become responsible for all inter-cell commu- 
nication. 

Before going into the details of the Riemann solver, we seek a simpler rep- 



resentation of (40 1. We start by defining the first term on the right hand side 



of ( 40 1 as the interior integral 

N^'^:=^fy^'HO.n<iHt,0)dt (41) 

and after rescaling, we define the last two terms as 

^i?+i/2 ■= ^^'He = +l).fHq\t,x,+y,)), (42a) 

<li/2 ■■= ^^"Hi = ~i)fHq\t,x,_,/,)). (42b) 

With these definitions in place, the discrete weak formulation can now be com- 
pactly written as a large MOL system: 



(fc) r^\ F(fc) 



j^Q['\t)= ivf(t)^ ^^[f;:/^,,,(t)-Frj-i/2(0j; (43) 

Interior Edge 

The only remaining piece to describe is how we compute the inter-cell flux 
values, f^{q''{t,x,_i/2)). 

4.I.I. The discontinuous Galerkin (DG) method: a choice of Riemann solvers 
Given that the representation of the solution, q^ , is not forced to be con- 
tinuous at the cell interfaces, care must be taken when defining the flux val- 
ues f^ {q'''{t,Xi_i/2)) ■ One could work with the so-called generalized Riemann 
solvers which take into account spatially varying function to the left and right 
of the discontinuity. One could also work with exact solutions for a classical 
Riemann problem in which one assumes constant states to the left and right 
of a discontinuity. In this work, we use Rusanov's method [58], which is an 
approximate Riemann solver that is commonly called the local Lax-Friedrichs 
(LLF) Riemann solver given its similarity to the Lax-Friedrichs method. Much 
like the Lax-Friedrichs method, this solver approximates each Riemann problem 
with two waves with equal speeds traveling in opposite directions, but LLF uses 
a local speed in place of a global speed to do so. 
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Evaluating the basis functions to the left and right of a single interface 
located at 2^^-1/2 provides artificially constant states, defined by a finite sum, 

Qr = Q['^- \/3Qf ^ + \/5Qf ^ - V7Q[*^ + ■■■ (44a) 

Qi - 0!i\ + VS Qf\ + ^/l Qf\ + V7 q\% + ■■■ . (44b) 

The LLF solver defines a single value based on these two constant states: 

/^ {Ql,Qr) ■■= I [(/ iQl) + f {Qr)) -a{Qr- Ql)] • (45) 

We use a local value of a, defined by a = max { | s^ | , | s^ | } , where s^ and s^ are 
the HLL(E) speeds [59] defined by 

s^ =niin|minA(^) (q\ , minA^P) iQi)\ , (46) 

s^ =max jmax A^P) (q) , maxA^^^ {Qr)\ , (47) 

and A'-P^ are the eigenvalues of /'. For this work, we take the simple arithmetic 
average Q — ^ {Qi +Qr), but we point out that Roe averages could also be 
used [53] , 

This completes the DG method of lines (MOL) discretization for our PDE. 
The only remaining part is to evolve the discrete coefficients through time. This 
is usually performed by explicit, high-order Runge-Kutta methods, resulting in 
a Runge-Kutta discontinuous Galerkin (RKDG) method. In principle, one may 
potentially work with this discretization to form a multiderivative integrator, 



which would require taking a second derivative of (43). However, difficulty will 



ensue when trying to compute the Jacobian of the complex system of ODEs in 



(40). Instead, we turn towards the Lax-Wendroff/Cauchy-Kowalewski discon- 



tinuous Galerkin methods to define a discrete second derivative. 

4-. 8. The discontinuous Galerkin (DG) method: a multiderivative formulation 

In this section, we extend the work of Qiu, Dumbser and Shu [S] to accom- 
modate multiderivative technology. An investigation into other fluxes would 
be an interesting topic of future study, including a comparison of various ap- 
proximate Riemann solvers using multiderivative technology [B]; moreover, an 
investigation into generalized Riemann solvers [HI HOI EI] may yield some inter- 
esting results. At present, our current goal is to lay the foundation that would 
be necessary for such investigations. 

Starting with the aim of defining a single stage value (or full update) for 



a multiderivative integrator, a DG implementation of equation (12) from S2.3 
requires the definition of q^tt as well as the definition of q\f. Here, we use 



q,tt = i.f'{q)fx)a; in place of q^tt = - {f'{<l)qt)x which was used in ^3.2 for the 
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second derivative. After factoring out a single spatial derivative, we can write 



([12|) as 

g - g" - AtLr + a*r - ^tpf'{q")f: - At(3*nq*)f:] , (48) 

where we have made use of the shorthand definitions 

/":=/(?"), r--^fiql- (49) 

We remark that after defining the appropriate spatial discretization, equation 



( 48 ) will look strikingly similar to the semi-discrete weak formulation already 



presented in equation ( 43 ) 



Before proceeding onward, we argue that a complete understanding of how 



to discretize ( 48 ) will provide the necessary building block for arbitrary multi- 



derivative Runge-Kutta integrators. For example, setting a* = /3* = 0, 



view this as a single-step method, and therefore equation (48) is nothing other 
than the Lax-Wendroff DG scheme already presented in the literature, provided 
the correct coefficients are inserted. Moreover, setting (5 — f3* — 0, this produces 
a two-stage Runge-Kutta method. Further stages can be included by adding in 
additional terms which introduces cumbersome notation. For clarity, we restrict 



our attention towards defining the discrete formulation of ( 48 ) which is general 
enough to accommodate more complicated multiderivative time integrators. 

The first step is to construct a Galerkin representation of the flux func- 
tion /'* as well as necessary coefficients for the second time derivative, g'^ := 
f'{q'^)f{q^)x- These two functions are the first and second, respectively, time 
derivatives of q before taking the final spatial derivative. The Galerkin coeffi- 
cients for these spatial derivatives can be constructed by taking advantage of 
the basis functions; 

Procedure 1.1 

INPUT: 

j Q,- > — a fist of Galerkin coefficients of q'^ . 
OUTPUT: 

a list of Galerkin coefiicients of f{q^), and 



{fI-'} 
{ 



Gl \ — a list of Galerkin coefficients of f'{q^)fx- 

1. For each element %, evaluate the flux function, fi{£) := f{q^{C)) ^^d the 
Jacobian Ji(^) :— /'(^f (C)) ^-t a list of M quadrature points, {^i, . . ■S,m}- 
Here, ^ is the canonical variable for grid element 7^ related to x through 
dSTl). 
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2. Compute the projection of / onto the basis functions using point values 
from Step 1: 

Ft^:={f,iO.'p'^''HO)- (50) 

This projection step produces a Galerkin expansion of the flux function 
on grid element %, that when written in the canonical variable is given 

by, 

M 

/.(0:=E^^'^W^^'^(0. (51) 

k=l 



3. Differentiate the flux function (51) on the interior of cell Ti to produce 



^^/^(0-£:E^''^| iO- (52) 



fc=l 



Note that this differentiation step is only valid on the interior of each cell. 
Evaluate this function at the interior quadrature points, {^1,^2, • ■ • ,^a/}- 
4. Define the Galerkin coefficients of g^ as, 

Gf :=(j.(e)-9./,(e),^(''^)(0)- (53) 

The required point values for the integration for Ji{S,m) were computed in 
Step 1, and the required point values of OxfiiCm) were computed in Step 
3 by differentiating the basis functions. 

We remark that differentiating the basis functions loses a single order of 
accuracy, but given that those terms involving g'' are multiplied by At = 0{Ax), 
we recover the desired order of accuracy. We are now prepared to describe the 
full multistage multiderivative DG method. 

Multistage multiderivative DG procedure 

1. Given a pair of Galerkin expansions g'* and q* , we compute four Galerkin 
expansions f^,g'^,f* and g* using Procedure 1.1. 

2. Define a modified flux function f^ via 

f'^:=af'^ + a*r+At{(3g'^ + /3*g*). 

We remark that we have a full Galerkin expansion of this modified flux 
function, that when restricted to a single element, is given by 



/''(-)^ = E^^'V('=no, 
fc=i 



20 



We can now simplify equation ( 48 1 by compactly writing it as 

q-g"-At/\. (54) 



3. Construct the time integrated weak formulation by multiplying (54) by a 
test function if^^' and integrate by parts over a grid cell 7^: 



■Ak) 



At r~ 



pC^) /-in\ T^W 



The only remaining piece is to define the flux values, as well as define 
a proper left and right fiux value for the Riemann solver. The interior 
integrals N^ are given by integration of the Galerkin expansion of /'': 

Once we have the Galerkin expansion of /'*, these integrals can be evalu- 
ated exactly. 

4. Solve Riemann problems for the modified flux function. In place of using 



equations (44a) and (44b) to define the left and right hand side Fi and 
Fr of the Riemann problem, we insert the extra time derivatives drawn 
from /'' and tuck them into this evaluation. That is, at the grid interface 
located at Xj_i/2, we redefine the left and right values of /^ to be the left 
hand side and right hand side evaluations of /'': 

M M 

F. :-/n-.^r/,) = T.^'-W''(^ = +1) = E VW^lFtl 

k=l k=l 

M M 

fc=i fc=i 

This completes the multiderivative description of the method within the 
discontinuous Galerkin framework. Extensions to multiderivative integrators 
that require more stages or more derivatives are a tedious, yet straight-forward 
extension of what has already been presented. Extensions to methods with more 
stages would require adding more values of a* and (3* to the definition of /''. 
Extensions to methods with more derivatives require bootstrapping previous 
Lax-Wendroff DG work ^, much like what has already been performed here. 
This would involve the appropriate extension of Procedure 1.1 using the Cauchy- 
Kowalewski procedure from equation (|4|, which would define an appropriate 
modification / of the flux function / that contains extra time information of 
the PDE. 

4-3. The discontinuous Galerkin (DG) methods: a choice of limiters 

The one detail that has been left out of this discussion has been the choice 
of limiting options. It is a well known fact that high-order linear methods 
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exhibit oscillatory behavior near discontinuities, and therefore to obtain phys- 
ically relevant results, one needs to choose a limiter that ideally retains high 
order accuracy in smooth regions, and reduces to first order accuracy locally 
at the location of the shock. In this work, we use the moment based limiter 
developed by Krivodonova [52] for our numerical simulations, although other 
limiters may certainly be used as well, with no change to the general framework 
presented thus far. This limiter is applied after each stage of the (multideriva- 

tive) Runge-Kutta method, and it modifies the higher order terms <Ql , fc > 2 ^ 
after each time step. One advantage of using multiderivative methods is that 
expensive limiters need to be applied less often when compared to high-order, 
single-derivative counterparts because of the reduction in the number of stages. 

5. Numerical examples 

In this section, we present results of the proposed method on a variety of hy- 
perbolic conservation laws. In §j ]5.1| - [5?2l we consider scalar examples: constant 
coefficient advection, and the Buckley-Leverett two-phase flow model. In §§5.3| 



- 5.4 we present results for the shallow water and Euler equations. All discon- 
tinuous Galerkin solutions were run using the open-source software DoGPack 
[63|. 

Unless otherwise noted, we use the following parameters: 

• The time integrator is the unique two-stage, two-derivative, fourth-order 



Runge-Kutta method (TDRK4) from equation (11). 



• All finite difference WENO simulations use the fifth order WEN05-Z re- 
construction from 93 and a constant CFL number oiv = 0.4 = max^* If'il*)] 

• Every DC result is fourth-order accurate in space, and uses a desired CFL 
number of j/ = 0.08. If the maximum allowable CFL number defined by 
I'max — 0.085 is violated, a smaller time step is chosen. For visualization, 
we plot exactly 4 uniformly spaced points per grid cell. 

For the two problems where other time integrators are compared against the 
TDRK4 method, we use coefficients in Table [T] for the DG simulations. 

5.1. Linear advection 

Our first example is a scalar linear advection equation with periodic bound- 
ary conditions: 

g,t+<?,:r=0, a;e[-l,l]. (55) 

5.1.1. Linear advection: a smooth example 

For a smooth example, we use initial conditions 

q{0,x) = qo{x) — sin(7ra;), (56) 
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At 



DG time stepping parameters 


V 


^inax 


Third-order (SSP-RK3) 


0.125 


0.130 


Low-storage SSP Runge-Kutta (SSP-RK4) 


0.44 


0.45 


Two-derivative Runge-Kutta (TDRK4) 


0.08 


0.085 



Table 1: CFL parameters used for DG simulations that compare time integrators. The SSP- 
RK3 method is the optimal third order SSP method developed by Shu and Osher [44) and 
described by Gottlieb and Shu |64| . SSP-RK4 is a fourth-order, low-storage method with ten 
stages developed by Ketcheson [39| . The maximum allowable CFL number, i-'max is near the 
maximum possible stable time step that each method permits for fourth-order spatial accuracy. 
We note that SSP-RK4 has a much higher CFL limit when compared to either TDRK4 or SSP- 
RK3 because it incorporates many more stages, but each stage requires expensive appli cati ons 
of limiters. It would be interesting to investigate optimized versions of TDRK4 using ( |12[ l as 
a building block that allow taking larger time steps without adding extra storage. 



and we run the simulation up to i = 2.0, at which point the exact solution is 
given by the initial conditions. This is one of two problems where we compare 
popular time integrators against the new method. Convergence studies are 
presented in Tables [2] and |3] For the WENO scheme and for this problem only, 
we run this problem with a large CFL number oi v — 0.9 in order to make the 
temporal error the preponderant part of the total error; accordingly, in Table 1 
we observe the fourth-order time accuracy of RK4 and TDRK4. 

Relative errors for the WENO scheme are defined by an L^ norm based on 
point-wise values: 



Error := 



AxE:"i(gr-g(t",:^.)r 



(57) 



Errors for the DG simulations likewise use a relative i^-norm, but this can be 
based on integration against the exact solution using moments of the solution 

(c.f. m)- 

5.1.2. Linear advection: a challenging example 

The second example we test is more challenging. The initial conditions are 
a combination of four bell-like shapes with different degrees of smoothness, that 
was originally proposed by Jiang and Shu |46) : 

r 1 



qo(x) 



[G{x, l3,z-6) + G{x, I3,z + S)+ 4G(x, /3, z)] , 

1, 

1- |10(a;-- 0.1)1, 

g [F{x, a,a — S) + F{x, a,a + 6) + 4:F{x, a, a)] , 

0, 



G(x,/3,z)=e-'3("-^)'; 



F{x, a, a) = ■\/max(l — a'^{x — a)^, 0). 



-0.8 <x< -0.6; 
-0.4 <x < -0.2; 
0< a; < 0.2; 
0.4 < x < 0.6; 
otherwise. 

(58a) 

(58b) 

(58c) 
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Mesh 


SSP-RK3 


Order 


RK4 


Order 


TDRK4 


Order 


25 


3.09 X 10-"=* 


— 


2.00 X 10-"4 


— 


1.52 X 10-""* 


— 


50 


3.79 X 10-°-* 


3.03 


9.68 X 10"°^ 


4.37 


8.77 X 10-"*' 


4.11 


100 


4.74 X 10-"^ 


3.00 


5.54 X lO-'J'^ 


4.13 


5.39 X lG-"7 


4.02 


200 


5.92 X 10-°6 


3.00 


3.37 X 10-°** 


4.04 


3.35 X 10-"8 


4.01 


400 


7.39 X 10""^ 


3.00 


2.09 X 10-"" 


4.01 


2.09 X 10-'"* 


4.00 


800 


9.24 X 10-°8 


3.00 


1.31 X 10-1" 


4.00 


1.31 X 10-1° 


4.00 


1600 


1.16 X lO-*^* 


3.00 


8.33 X 10-1^ 


3.97 


8.33 X 10-1" 


3.97 



Table 2; Advoction equation: smooth example. Convergence analysis for the WENO scheme 
applied to the advection equation \55\ with periodic boundary conditions and smooth initial 
conditions ( |56| . The table shows a comparison of the relative L^-norm of the errors in the 
numerical solutions obtained with WEN05-Z spatial discretization and different time integra- 
tors: 3rd-order strong stability preserving Runge-Kutta (SSP-RK3) |44| . classical 4th-order 
Runge-Kutta (RK4), and the new two-derivative 4th-order Runge-Kutta (TDRK4) method. 
The Courant parameter chosen for this problem was a constant CFL oi u = 0.9 = At/ Ax. 
In order to observe fourth-order accuracy for the fourth-order methods, we needed to in- 
crease the CFL number and therefore increase the temporal error. Results for smaller CFL 
numbers reached machine precision before announcing their accuracy, where both RK4 and 
TDRK4, indicated convergence orders between 4*'' and 5"" order. The 'Order' columns refer 
to the algebraic order of convergence, computed as the base-2 logarithm of the ratio of two 
successive error norms. We remark that for all CFL numbers and weighting schemes tested, 
including WENO-Z, WENO-JS and linear weights, both time integrators RK4 and TDRK4 
have comparable errors. 



Mesh 


SSP-RK3 


Order 


SSP-RK4 


Order 


TDRK4 


Order 


5 


1.17 X 10-°^ 


— 


6.18 X 10-°'! 


— 


8.21 X 10-°* 


— 


10 


1.32 X 10-°i 


3.16 


3.87 X 10-°^^ 


4.00 


5.99 X 10"°^ 


3.78 


20 


1.60 X 10-°^ 


3.04 


2.43 X 10-°^ 


3.99 


3.91 X 10-°'' 


3.94 


40 


1.99 X 10-°*^ 


3.01 


1.52 X 10-°^ 


4.00 


2.47 X 10-°^ 


3.98 


80 


2.48 X 10-°^ 


3.00 


9.51 X 10-°''' 


4.00 


1.55 X 10-°** 


4.00 


160 


3.10 X 10-°« 


3.00 


5.94 X 10-1° 


4.00 


9.69 X 10-1° 


4.00 


320 


3.87 X 10-°** 


3.00 


3.72 X 10-11 


4.00 


6.03 X 10-11 


4.01 


640 


4.84 X 10-1° 


3.00 


2.24 X 10-1" 


4.05 


4.39 X 10-1" 


3.78 



Table 3: Advection equation: smooth example. We present DG-results comparing the new 
multiderivative scheme TDRK4 against third-order SSP-RK3 [441 164) and a state of the art 
low-storage fourth-order Runge-Kutta method (SSP-RK4) I39| . The CFL numbers chosen for 
each scheme are presented in Table l!] which are near the maximum stable CFL limit for each 
scheme. 
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The constants are a = 0.5, z = -0.7, S = 0.005, a = 10 and /3 = logio(2)/(36 S^). 
With this choice of constants, the Gaussian curve centered at a; = —0.7 contains 
two discontinuities, and the bump centered at x — 0.5 also has two discontinu- 
ities. Moreover, this last bump has no derivative at a; = 0.405 and x = 0.595. 

After testing many different reconstruction procedures of varying orders, 
many different time integrators and a large range of parameters, our observations 
indicate that the WENO methods are sensitive to the parameters used for this 
test problem. It has already been observed that WENO schemes may require 
tuning in order to achieve good results [5T]; this is not the subject of this work, 
but we find it necessary to point out that WENO simulations do not emulate 
robust behavior for this problem. A numerical investigation into their behavior 



on linear problems with difficult initial conditions such as (58a)- (58c) would be 
interesting. 

Results are presented in Figure [T] The under-resolved DG simulations tend 
to be quite diffusive when compared to the WENO schemes, but behave much 
more predictably. 

5.2. Buckley-Leverett 

The Buckley-Leverett equation is a non-linear scalar problem with a non- 
convex flux function. Years past its invention, it has become a standard bench- 
mark problem for many hyperbolic solvers. The model describes two-phase flow 
through porous media, where the application is oil recovery from a reservoir 
containing an oil- water or oil-gas mixture of fluids f65] . In rescaled units, the 
flux function for this problem is defined through a single free parameter M with 

We use AI = 1/3, and we take the computational domain to be [—1,1] with 
initial conditions prescribed through 

.(0,.) = ^' if-V2<-<0, (,,) 

I 0, otherwise. 

For small time values, the exact solution to this problem can be found by 
solving two Riemann problems, where each pair of states produces a typical 
'compound wave' that is the combination of a leading shock and a trailing 
rarefaction. Each shock propagates with a speed dictated by the Rankine- 
Hugoniot conditions 

/'(,■) ^ Hq^lM, 

where Qc is the constant value that does not change. For the Riemann problem 
located at x = —0.5, we have Qc = 1, and q* ~ 0.1339745962155613, and 
for the Riemann problem located at a; = 0.0, we have qc — and q* = 0.5. 
Characteristics between q* and qc propagate with speed f'{q), and therefore 
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Figure 1: Advection equation: a discontinuous example. Results for the two-derivative 
(TDRK4) method for WENO and DG spatial discretizations. Simulations are run to a fi- 
nal time of i = 8.0 for a total of four full revolutions. The two images on the right are zoomed 
in images of the top of the square wave. WENO simulations use a CFL number oi u = 0.4, 
and the DG simulations use a constant CFL number of i^ = 0.08. We again plot four points 
per cell for the DG solutions. 
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Figure 2: Buckley-Leverett double Riemann problem: WENO solutions. WENO results at 
the final time of t = 0.4. All time integrators compared for this method are giving similar 
results. We remark that for this problem, the WENO schemes tend to be more diffusive where 
the rarefactions form when compared to the DG schemes. 



to fill the rarefaction fan, we plot a range of values (s + tf'{q),q), where g is a 
sampling of values between Qc and q* , and s G {—0.5, 0.0} is the location of each 
Riemann problem. Results for the two spatial discretizations are presented in 
Figures [2] and [3] at a final time of i = 0.4. 

As a final note, we remark that for this problem only, we modify our scheme 
to deal with some pathological issues. Given that the flux function is non- 
convex, the WENO simulations use the analytical value for a, which is approx- 
imately a := maXqg[o^i] \f'{q)\ ~ 2.205737062. For the DG simulations, we use 
the HLL(E) Riemann solver [59] . In this case, the HLL(E) solver performs bet- 
ter than the local Lax-Friedrichs (LLF) solver given that f'{q) > for all q in 
our domain. The LLF solver approximates this problem with two waves trav- 
eling in opposite directions, whereas the exact solution travels in one direction 
only, which the HLL(E) correctly accounts for. 
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Figure 3: Buckley-Leverett double Riemann problem: DG solutions. DG results at the final 
time of t = 0.4. We plot exactly four uniformly spaced grid points per cell. Tfie top row uses 
mx = 40 grid cells, and the bottom row uses rrix = 80 grid cells. The two frames on the right 
are zoomed in images of the solution near x = 0.0, one of the few places where we were able 
to find differences in the solutions. We note that the coarse solution is under resolved because 
it only has a single cell to capture the structure at the top of the solution. CFL numbers 
are chosen as in Table IT] which are close to the maximum stable time step allowed. All time 
integrators arc qualitatively giving the same result. 
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5.3. Shallow water 

The shallow water equations define a hyperbolic system with two equations 
for water height h, and velocity u. The hyperbolic conservation law that these 
two variables satisy is given by 



h \ I hu 

hu J I hu^ + \gh 



2 , 1.U2 I =0. (61) 



For our simulation, we take g = 1. We demonstrate our method on the dam 
break Riemann problem [66j with initial conditions defined by 

1(1,0) , otherwise. 

We use a computational domain of [0, 1] with outfiow boundary conditions, and 
stop the simulation at a final time of i = 0.2. Results for the two-derivative 
time integrators are presented in Figure|4]for the WENO method, and in Figure 
[5] for the discontinuous Galerkin scheme. Exact solutions for this problem can 
be found in textbooks [66 . 

5.4. Euler equations 

The Euler equations describe the evolution of density p, momentum pu and 
energy 8 of an ideal gas: 



(63) 



where p is the pressure. The energy £ is related to the primitive variables p, u 
and p by 

e P 1 2 

7 — 1 2 

where 7 is the ratio of specific heats. For all of our simulations, we take 7 = 1.4. 

5.4.1. Euler equations: a shock tube Riemann problem 

We present a classic test case of a difficult shock tube, which is commonly 
referred to as the Lax shock tube. The initial conditions are those defined by 
Harten [IZIEH]: 



p ' 




' pu 


pu 


+ 


pu^ +p 


£ 1 


■t 


V {^+P)u 



^T _] (0.445, 0.3111, 8.928)^, if x < 0.5, 
(0.5,0,1.4275)"^, otherwise. 



{p,pu,£) = S ,„ , „ , ,„^,,7. _^,, ,_ (64) 



Exact solutions for this problem are well understood, and there are many text- 
books that describe how to construct them [66l IM] • For this set of data, the 
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Figure 4: Shallow water Riemann problem: WENO solutions. Here we present two resolutions 
on top of each other using the new multiderivative scheme presented in j3.2| We use rrix = 50 
points for a coarse resolution and nix = 150 points for a finer solution. Top row: water height 
h at the final time with a zoom in of the rarefaction to the right. Bottom row: momentum 
hu, with a zoom in of the shock to the right. 
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Figure 5: Shallow water Riemann problem: DG solutions. Here we present two resolutions on 
top of each other using the new multiderivative scheme presented in j ]4.2[ We use mx = 20 
cells for a coarse resolution and m^ = 60 cells for a finer solution, and we plot exactly four 
uniformly spaced points per cell. Top row: water height h at the final time with a zoom in of 
the rarefaction to the right. Bottom row: momentum hu, with a zoom in of the shock to the 
right. 
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solution contains a left rarefaction, a contact discontinuity and a shock wave 
traveling to the right. 

We select t = 0.16 for the final time of our simulation [5], and we use a 
computational domain of [0, 1] with outflow boundary conditions. Results for 
this problem are presented in Figures [6] and [7] For the WENO simulations, we 
additionally compare the two-derivative two- derivative time integrator TDRK4 
against the SSP-RK3 with two different CFL numbers, ly — 0.01 and then 
v = 0.4. We present results for the time integrator comparison in Figure [Sj 

5.4-2. Euler equations: shock entropy 

Our final test case is another problem that is popular in the literature [15] . 
The initial conditions are 

(p, w, p) = (3.857143, 2.629369, 10.3333) , x < -4, 

(p,-u,p) = (1 + esin(5a;),0,l) , x>-A, 

with a computational domain of [—5,5]. The final time for this simulation is 
t — 1.8. With e = 0, this is a pure Mach 3 shock moving to the right. We follow 
the common practice of setting e — 0.2. 

Results for WENO simulations are presented in Figure [9J and DG results are 



presented in Figure 10 For a reference solution, we plot a WENO simulation 
that uses the SSP-RK3 method described in Gottlieb and Shu ^4], with ttt.^ 
6000 points and a small CFL number of :^ = 0.1. 



6. Conclusions and future work 

We have presented an alternative to popular high-order time integration 
methods for hyperbolic conservation laws. The explicit multistage multideriva- 
tive time integrators that are the subject of this work sit halfway in between 
explicit (single-derivative) Runge-Kutta methods and high-order (single-stage) 
Taylor methods. Given the physical locality of hyperbolic problems, multi- 
derivative integrators can be implemented with a low-memory footprint and 
low message passing overhead, which makes them attractive in light of mod- 
ern day computer architecture. Numerous numerical examples were presented 
that included multiderivative schemes for high-order discontinuous Galerkin and 
WENO methods. In order to implement the multiderivative technology, we 
leveraged recent work on Lax-Wendroff type time integrators for the two spatial 
discretizations investigated here. Numerical results for the new multiderivative 
schemes are promising: they are demonstrably comparable to those obtained 
from popular high-order SSP integrators, yet given their portability, locality 
and naturally low- memory footprint, we may see them gain popularity in the 
near future. 

Future work will focus on exte nsions to higher dimensions, as well as an 
investigation into numerical properties of multiderivative schemes. We would 
like to explore developing multiderivative methods with SSP properties, as well 
as low-storage 'many'-stage variations of the two-derivative method presented 
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Figure 6: Shock-tube Riemann problem; WENO solutions. Shown here are WENO simula- 
tions with Harten's initial conditions ( |64[ l for the shock tube problem. We plot observable 
quantities from top to bottom: density p, velocity u, and pressure p. Left columns are the full 
solution, and right columns are zoomed in parts of the same data points. We present a coarse 
solution with m^ = 100 points and a finer solution of rrix = 300 points. The zoomed in image 
of the density indicates the top part of the square section between the contact discontinuity 
and the right traveling shock. The zoomed in images for the velocity and pressure focus in on 
the right foot of the rarefaction fan. 
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Figure 7: Shock-tube Riemann problem: DG solutions. We present two DG simulations for 
the physical observables, which from top to bottom are: density p, velocity n, and pressure 
p. Left columns are the full solution, and right columns are zoomed in parts of the same data 
points. We present a coarse solution with rrxx = 70 grid cells and a finer solution of rux = 210 
grid cells, and we plot four uniformly spaced points per cell. The axes are identical to those 
in Figure |6] 
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Figure 8: Shock-tube Riemann problem: WENO solutions. Shown here are WENO simu- 
lations comparing two different time integrators: the third-order SSP method of Osher and 
Shu against the new two-stage, two-derivative method (TDRK4). We repeat that all WENO 
simulations with the TDRK4 time integrator use a CFL number of f = 0.4. Here, we observe 
that the third-order SSP method with CFL numbers of i^ = 0.4 and v = 0.01 behaves quali- 
tatively the same as the TDRK4 method. We plot observable quantities from top to bottom: 
density p, velocity m, and pressure p. Left columns are the full solution, and right columns 
are zoomed in parts of the same data points. The resolution for this picture is the coarse 
solution of rria; = 100 identical to the coarse solution used in Figure |6] Again, we remark that 
the accepted time integrator is in close agreement with the proposed method. 
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Figure 9: Shock-entropy: WENO solutions. Shown here are WENO solutions to the shock- 
entropy interaction problem presented in j ]5.4.2| We plot observable quantities from top to 
bottom: density p, velocity n, and pressure p. Left columns are the full solution, and right 
columns are zoomed in parts of the same data points. We present a coarse solution with 
rrta; = 100 points and a finer solution of m^ = 300 points. 
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Figure 10: Shock entropy. Shown here are DG solutions to the shock-entropy interaction 
problem presented in §5.4.2[ We use a coarse mesh of rrix = 70 and a finer mesh of rux = 210. 
The reference solution is identical to the one used in Figure |9] that is a WENO simulation 
with 6000 points, and CFL number of i/ = 0.1 with the third-order SSP-RK3 time integrator. 
Again, the DG schemes tend to be more diffusive on this problem when compared to the 
WENO method, and the limiter used in this work tends to clip the peaks for the coarse 
resolutions. We repeat that we are plotting exactly four points per cell, and therefore, the 
observable oscillations for the coarse resolution are happening on a sub-cell level. Although 
not shown, plots of cell averages only do not visibly demonstrate these oscillations. Results 
for other time integrators are found to be comparable. 
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in this work. Additionally, we would like to investigate the utility of multistage 
multiderivativc methods for solving parabolic partial differential equations. 
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